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AN ABSOLUTELY STABLE hp -HDG METHOD FOR THE 
TIME-HARMONIC MAXWELL EQUATIONS WITH HIGH WAVE 

NUMBER 

PEIPEI LU, HUANGXIN CHEN, AND WEIFENG QIU 


Abstract. We present and analyze a hybridizable discontinuous Galerkin (HDG) method 
for the time-harmonic Maxwell equations. The divergence-free condition is enforced on 
the electric field, then a Lagrange multiplier is introduced, and the problem becomes the 
solution of a mixed curl-curl formulation of the Maxwell’s problem. The method is shown 
to be an absolutely stable HDG method for the indefinite time-harmonic Maxwell equations 
with high wave number. By exploiting the duality argument, the dependence of convergence 
of the HDG method on the wave number k, the mesh size h and the polynomial order p is 
obtained. Numerical results are given to verify the theoretical analysis. 


1. Introduction 


The time-harmonic Maxwell boundary value problem reads as follows: 

curl curl u — k 2 u = f in G, 

curl u x n + i ku 1 = g on dfl, 


(l.la) 

(1.1b) 


where C M 3 is a bounded, uniformly star-shaped polyhedral domain, the wave number 
k is real and positive, i denotes the imaginary unit, n denotes the unit outward normal 
to dfl, and u l = (n x u) x n denotes the tangential component of the electric field u. 


Equation (1.1b) is the standard impedance boundary condition which requires g ■ n — 0, 


thus, g = g. The above Maxwell equations are of considerable importance in the engineering 
and scientific computation. In this paper we assume the current density is divergence-free 
(namely div/ = 0), hence the electric field u is also divergence-free. 

The Maxwell’s operator is strongly indefinite for high wave number k, which brings difficul¬ 
ties both in theoretical analysis and numerical simulation. Various numerical methods which 
include finite element methods (FEM) [2011211 H2l EEIl [7J [26] , discontinuous Galerkin (DG) 
methods [23 [2H El uni uni [321 [in [10] and weak Galerkin FEM method [12] have been de¬ 
veloped to solve the Maxwell’s problem. In particular, Feng and Wu ra recently proposed 
and analyzed an interior penalty discontinuous Galerkin (IPDG) method for the problem 
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(1.1) with high wave number, which is uniquely solvable without any mesh constraint. DG 


methods have several attractive features which include the capabilities to handle complex 
geometries, to provide high-order accurate solutions, etc. But the dimension of approxima¬ 
tion DG space is much larger than the dimension of the corresponding conforming space. 
Hybridizable discontinuous Galerkin (HDG) methods [5j were recently introduced to ad¬ 
dress this issue. The HDG methods retain the advantages of standard DG methods, and the 
resulting system is only due to the unknowns on the skeleton of the mesh. 

Two HDG methods were presented in [22 J for the numerical solution of the Maxwell 
problem. The first HDG method enforces the divergence-free condition on the electric field 
and introduces a Lagrange multiplier. It produces a linear system for the degrees of freedom 
(DOF) of the approximate traces of both the tangential component of the vector field and 
the Lagrange multiplier. The second HDG method does not enforce the divergence-free 
condition and results in a linear system only for the DOF of the approximate trace of the 
tangential component of the vector field. Compared to the IPDG method for the time- 
harmonic Maxwell equations in [16, [TO] , the two HDG methods have less globally coupled 
unknowns. The well-posedness, conservativity and consistence of the two HDG methods, 
together with a numerical demonstration, have been shown in [22J. However, no convergence 
analysis is given in [22]. Recently, the /?,-convergence analysis of the second HDG method 
was considered in HU- In this paper we are interested in the /ip-convergence analysis for 
the first HDG method mentioned in ra which solves a mixed curl-curl formulation of the 
time-harmonic Maxwell equation 


curl curl u — k 2 u + Vd = f 
div u = 0 
curl u x n + inu* = g 

a = 0 


in ff, 

(1.2a) 

in D, 

(1.2b) 

cf 

Pl 

o 

(1.2c) 

on (Kl. 

(1.2d) 


where a is a scalar Lagrange multiplier used to enforce the divergence-free condition. Taking 


the divergence of the equation (1.2a) yields Act = 0, which together with the boundary 


condition (1.2d) implies that a = 0 throughout the domain. Hence, under the divergence- 


free condition of the current density, the equations (1.1) and (1.2) are equivalent. 


We aim to develop an HDG method which is absolutely stable without any mesh constraint 


for the above mixed curl-curl formulation (1.2) and reveal the dependence of convergence for 


the HDG method on the wave number k, the mesh size h and the polynomial order p. We 
mention that only simple L 2 -projections are used in our analysis which is different from the 
projection-based error analysis in [6], and the p-dependence of the stability estimate and the 
convergence can be derived. We also mention that the stabilization parameters in our HDG 
method are different from those in ra- The focus of our analysis is to apply the duality 
argument to establish the rigorous stability estimate and error analysis for the HDG method 


proposed for the mixed curl-curl formulation (1.2). Intrinsically, the regularity estimate of 


the solution of the dual problem used in this paper can be obtained due to introduction of 
a Lagrange multiplier in the mixed curl-curl formulation. This is also the reason why the 
Helmholtz decomposition technique can be avoided in the analysis and the p-estimate can 
be derived. We first apply the duality argument to obtain the estimates for || || ^ and the 

error ||u — Uh]\j h , then the estimates for other variables of the HDG method can be further 
obtained. Up to our best knowledge, we give the first p-estimate of numerical methods using 
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piecewise polynomial solution spaces for solving the time-harmonic Maxwell equations with 
high wave number. 

The remainder of this paper is the following. We give some notations, introduce the HDG 
method for the mixed curl-curl formulation of the time-harmonic Maxwell equations (1.2) 
and present the main results of stability estimates and error estimates in the next section. 
Section [3] and section [| are devoted to providing detailed proofs of the stability estimates 
and error estimates respectively. In section [5j we discuss the stability estimates and error 
estimates for the HDG method under some ideal assumptions of the problem (1.1) and the 
associated dual problem. In the final section, we give some numerical results to confirm our 
theoretical analysis. 


2. Notation, HDG method and main results 


Let f := —if, a := —id and g := —i g. The HDG scheme for the equation (1.2) is based 
on a first-order system of this equation, which can be rewritten in a mixed formulation as 
follows: 


iw — curl 14 = 0 in D, ( 2 . 1 a) 

curl w + i K 2 u + Va = f in D, ( 2 . 1 b) 

divu = 0 in D, ( 2 . 1 c) 

w x n + ku* = g on <9h2, ( 2 .Id) 

cr = 0 on dQ. ( 2 -le) 


Throughout the paper we use the standard notations and definitions for Sobolev spaces 
(see, e.g., Adams CD- We denote by 7h a conforming triangulation of D made of shape- 
regular simplicial elements. We denote by h F the diameter of T G 7h and h = max^er^ hr, 
the collection of faces is denoted by 8 ^, with the collection of interior faces by £° and the 
collection of boundary faces by £f, the collection of element boundaries by d7 h : = {<9T|T G 
T/ t }. We let C denote a positive number independent of the mesh size, polynomial order 
and wave number, whose value can take on different values in different occurrences. The 
corresponding finite element spaces for the HDG method for the first-order system (2.1) are 
defined to be 


V h := {r G L 2 (Q) : r\ T G P p (T),VT G 7 h }, 

U h := {u G £ 2 (D) : v\ T G P p (T),VT G 7 h }, 

Ml := {r 7 G L 2 (E h ) : 77 ^ G P p (F), (77 • n)\ F = 0, VP G £4, 

Q h := {q G L 2 (D) : q\ T G P P (T),VT G 7 h }, 

M h := G L 2 (£,) : £\ F G P p (P),VP G £*}, 

where the polynomial order p > 1, X 2 (D) = [L 2 (fl)] 3 , X 2 (£h) = [L 2 (£/ l )] 3 , P p (T) = [P p (T )] 3 
and P p (P) = [P p (F)] 3 . Here, P p (D) denotes the space of complex-valued polynomials of 
degree at most p on D. Let Pm denote the standard L 2 -projection operator from L 2 (£/ l ) 
onto P p (£/ l ). In addition, we set M h (g ) := {£ G M h : £ = P M g on <9fl}. Similarly, P M 
denotes the standard L 2 -projection operator from T 2 (£/ l ) onto P p (£/ t ). We use IIv- II u, 
to denote the standard L 2 -projection onto Vh,Uh and Qh respectively. I 11 the analysis, we 
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shall use the following approximation results of L 2 -projections: 


\w 


it 


u v w\\ 7h < civ/p t \\w\\ t p 
Huu\\ Jh < Ch s /p s \\u\\ Si u 
o-Yl Q a\\ Jh < Ch^/p^\\a\\ p>n 
- IIvtw||o,5T < Ch l 2 /p t 2 ||tn|| ti T 

— Huu\\o,dr < Ch s 2 /p b 2 ||u|| S) t 

IIcr - riQ(T||o,9r < C/Z^/fa^IHUr 
w - P M w\\ 97h < C'/i t- 3/p* - 2||iw|| t) n 
||w - PmU llao-fc < Ch s ~^/p s ~^\\u\\ Sj Q 
Pmct Hao-fc < CTi^/P^IMkn 


\w 


u 


\cr 


Here 


0<f <p+l, (2.2a) 

0 < s < p + 1, (2.2b) 

0<P<P+1, (2.2c) 

WT e7 h , 0 < t < p + 1, (2.2d) 

VT e T h , 0 < s <p + 1, (2.2e) 

VTeTfc, 0<P<p+l, (2.2f) 

0<t<p+l, (2.2g) 

0 < s < p + 1, (2.2h) 

0 < /? < p + 1. (2.21) 

k = (Etct,, II' llo.r)* and II' liar* = (Ets^ II' llopr)^- The above results hold due to 
the hp approximation theory of polynomials and trace inequality when T/, consists of shape- 
regular simplices (cf. [25j [9, 0, 5) HU Gj). The above /r-dependence approximation results 
hold when 7h consists of shape-regular polyhedral elements. Thus when we only consider 
k- and /r-dependence in our analysis, T/ t can be a conforming mesh consisting of shape- 
regular polyhedral elements. This is due to the fact that only the approximation results 
(2.2a)-(2.2c) have been deduced recently in the literature (cf. [2]) when the mesh consists of 
general polyhedral elements. The p-dependence of convergence for the trace estimate of the 
polynomial L 2 -projection (cf. (2.2d)-(2.2f)) was first studied in |4] on simplicial element, and 
as far as we know, no extension of the estimates (2.2d)-(2.2f) to the L 2 -projection defined 
on the general polyhedral element has been obtained. 

We define the bilinear forms 

fa, Ck := fa, Or, fa, 0&i h ■= fa’ & 9T 

TS 7 h T&7 h 

fa, Ch h ■■= fa’ Or, (p, ()<n h ■= fa’ 

T&7 h T&7 h 


where ( 77 , C)r> (respectively, (p,C )d) denotes the integral of 77 ■ £ (respectively, p£) over 
Del 3 and (t],C)d (respectively, (p, C )d) denotes the integral of 77 (respectively, r]() over 


D C M 2 . 


The HDG method for the first-order system (2.1) yields a solution (wh,Uh,u h , <Th,&h) £ 
Vh x Uh x M* h x Q h x M h { 0) such that 


(i w h , r h ) 7h - ( u h , curl r h ) 7h + {u\ x n, r h ) d7fi = 0 (2.3a) 

(w h ,curlv h ) 7h - (wl x n,v h ) d7h + (i/s - 2 u h ,v h ) 7h 

- (op, div v h ) 7h + (a h ,v h ■ n) d7h = (f,v h ) 7h , (2.3b) 

- ( u h , Vq h ) 7h + (ul ■ n, q h ) d7h = 0, (2.3c) 

x n, t]h)d7 h \an = 0, (2.3d) 

(w[ x n, n h )an + (««!,, »7 h )an = fa, n^asu (2.3e) 

W • n, t h ) a7h = 0, (2.3f) 
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for all (r h ,v h ,T] h ,q h ,^ h ) EV h xU h x xQ h x M h (Q), where 

w h = w h + T t (u t h -u t h )xn, u n h = ul + T n (a h -a h )n. (2.4) 

Here, for any vector r E M 3 , r n = (r • n)n denotes the normal component of the vector 
r. The parameters r t and r n are the so-called local stabilization parameters which have an 
important effect on both the stability of the solution and the accuracy of the HDG scheme. 
We choose r t = p/h and r n = (1 + n)h/p in this paper. 


Remark 2.1. The mixed curl-curl formulation (1.2) can also be applied to the Maxwell 
equations (1.1) with div / ^ 0. In this case div U = 6 ^ 0 with 0 a given variable. Indeed, 
taking the divergence of the equation ( 1.1a) im plies that 6 satisfies that — n 2 6 = div/. 
Then taking the divergence of the equation (1.2a) again yields Ad = div / + n 2 6 = 0, which 
together with the boundary condition (1.2d) also implies tha t d = 0. Hence, the HDG scheme 
in this paper can also be used for the Maxwell equations (1.1) with div / ^ 0. However, 
we mention that if the HDG method is used with non divergence-free current density, the 
regularity estimates in [13] can not be applied. Thus, the theoretical analysis throughout 


this paper holds only under the assumption of divergence-free current density in (1.1a). 


When / is divergence-free and g E ii/(<9D), the solution of the first-order system (2.1) 
satisfies that u E H^ +a (Q) and w E i? 5+ "(y2) ; and there holds (cf. [13]) 

«IMI!+«,« + IMIi+ a ,n < C(l + K)M(f,g) + C|| 0 ||i i an . (2.5) 

To state our main results, we need a regularity assumption of the dual problem. Let and 
(p be the solution of the following dual problem: 


curl curl th — k 2 T r + Y\7p 
div \I> 

curl tL x n — 

V 


iJ 

in f2, 

(2.6a) 

0 

in 

(2.6b) 

0 

on dfl, 

(2.6c) 

0 

on dfl, 

(2.6d) 


where VJ E Uh C L 2 (Q). Due to the fact that is a bounded, uniformly star-shaped 
polyhedral domain, the solution (T r ,(^) has the following regularity estimate (cf. [iT3l [ID] ): 

«ll 1 ® r ll|+a I n + ll curl ^lli+a,n + + K )ll^ r llo,n + (1 + «)||V(^||o,o < C{ 1 + k)\\ J|| 0) n, (2.7) 

where 0 < a < a and 0 < a < ^ is a parameter only depending on When is convex, 
the above estimate holds true for all 0 < a < |. Moreover, when is also a C 2 star-shaped 
domain, (2.7) holds true for a = |. In the following, we show that (2.7) holds true under 
the assumption of the domain It is easy to see that th satisfies 

curl curl th — k 2 ^ = i(J — Vtp) in D, (2.8) 

div T' = 0 in D, (2.9) 

curlT' x n — = 0 on dfl. (2-10) 


By (2.8), we have 

(curl curl T' - Vq) n = (i(J - Vp), Wq)n 
By doing integration by parts and ( |2.9 ), we have 

(J-Vp,Vq) n = 0 V q E G/ 3 (f2). 


V q E C™(Q). 
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By density argument, we have 

(J ~ V(/p, Vg)n = 0 VqeH^tt). 

We easily obtain div (J — Vp) = 0 and || Vy?||o,o < || J||o,n- So, we can conclude that the 
estimate (2.7) holds true when hi is a bounded, uniformly star-shaped polyhedron (cf. [13]). 

Now we are ready to outline the main results in the following by showing the stability 
estimates of the discrete solutions from the HDG method (2.3) and the associated error 
estimates. 


Theorem 2.1. Let w , u. a) and (wh,Uh,u h , <Jh,&h) solve the equations (2.1) and (2.3). 
We assume that (2.1) holds with 0 < a < a < | and ( w,u ) € H^ +a (Q) x H^ +a (Q). Then 
the HDG method (2.3) is absolutely stable. When nh/p > Go, we have 


\9w 

'stab 11 

K 

l 


IWk<c(cS, b ll-llofl + c. 

< cU- + G s 2 tab )||/||o,Q + (— + Gstab) Iloilo, on), 

\ hi K2 / 

ll^hllaT h < G((— )5 +p/r _ 2) rG s 2 tab ||^||o,o + Gstabll^llo.ao 

T) \ rC rC 


P 

7 s~i -1 . Ac(l +K,)h 2 . K,(l-\-K,) 2 h 2a ~^ 1 . (1 +K) 2 h 2c * 

where C sta b := 1 + v p2 J + v p2 a+i -+ —• 


( 2 . 11 ) 

( 2 . 12 ) 

(2.13) 


Theorem 2.2. Let (w, u, a) and (wh,Uh,u h ,ah,&h) solve the equations (2.1) and (2.3). 
We assume that (2.1) holds with 0 < a < a < \ and ( w,u ) £ H^ +a (Cl) x i 3 " 2 +Q (n). When 
nh/p > G 0 , we have 


\u 


\w 


u h\\7 h 5s C{R W ||m|| i +Qj q + R u ||w|| i+q,^) > 


w h \\ 7 ,<C((()i+° + K.R 


in 


P 




+ C(nR u + k(1 + (1 + hi) 2 )(-) 2+a 

\ p 




(2.14) 


(2.15) 


where R n „ := (1+ ^I+i + + (1+h '- >2 \ —- and R. 


■= (l + (l + K)^)(^)^ + (l + «+(l + «)^)(J) 


h\2a 


Remark 2.2. For the solutions of the first-order system (2.1) which admit the regularity as 
in (2.5), when nh/p < Go, one may tune the parameters Tt and r n (cf. Remark 3.1) and also 
get the stability estimates and error estimates for the discrete solutions of the HDG method 
(2.3). When we consider only n- and ^.-dependence, the above results hold when “J h consists 
of general polyhedral elements. 


3. Stability estimate 


In this section we shall show that the HDG method (2.3) is absolutely stable. We first 
present a lemma which shall be used to estimate the stability estimate of u^. 

Lemma 3.1. Let (w^, Uh, u/, <Th., cp,.) be the solution of the problem (2.3). It holds that 

WhHK - u h )\\ 2 dJh + || rj(a h - a h )\\ 2 dJh + ^||uJo, aa < ||/|MMk + (3- 1 ) 


\\™h\\l h < ^\\ u h\\l h + 2||/|| 0) n||«h|k + -Iloilo,an- 

rv 


(3.2) 
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Proof. We first choose r h = w h ,v h = u h: t] h = u h ,q h = cr h ,£ h = &h in ( |2.3a| )-( |2.3f[ ) to get 
the following equalities: 

(i w h ,w h ) 7h - (u h ,curlw h ) 7h + (u[ x n,w h ) d7h = 0, 

(curl w h ,u h ) 7h + (w\ x n,u h ) d7h - x n,u h ) d7h 

+ (i n 2 u h , u h ) 7h - (a h , div u h ) 7h + (a h , u £ ■ n) d7h = (/, u h ) 7h , 

(div u h , a h ) 7h - « ■ n, of) d7h + (v% ■ n, a h ) d7h = 0, 

{™h x n ’ «h)e T h \osi = 0, 

(w h x n, u h ) dn + (KU h , u h )an = (g, u h )dn, 

(u n h ■ n,d h ) d7h = 0 , 


(3.3a) 

(3.3b) 

(3.3c) 

(3.3d) 

(3.3e) 

(3.3f) 


where (3.3b) and (3.3c) are obtained by integration by parts. Furthermore, noting the 
definitions of Wh in (2.4) and applying complex conjugation to (3.3a), (3.3c) and (3.3f), we 
get the following equalities after simple manipulations: 

- (i w h ,w h ) 7h - (curl w h ,u h ) 7h - ( T t (u{ - u h )\u h ) d7h + (ku^u^qq = {g,u h ) gn , 

( cr h , div u h ) 7h - {cr h , u n h • n) d7h + (a h , ■ n) d7h = 0, 

- ($h, K ■ ™)d7 h = 0. 


u h V = ui 


Adding the above three equalities and (3.3b) together and noting that (u l h 
we have 

- (i w h , w h ) 7h + (r t « - u h ), u{ - u h ) d7h + (KU h , u h ) dn 

+ (r n {a h - d h ), °h ~ Oh)d7 h + (iK 2 u h , u h ) 7h = (/, u h ) 7h + (g, u h ) dn , 
which implies the lemma by the Cauchy-Schwarz inequality. 


u 


h"> 


□ 


Next we shall utilize the dual argument to give the L 2 -norm es tima te of u^. Given 
u h G L 2 (Q), we introduce the first-order system of the dual problem (2.6) with J = u h : 


i<f> — curl th = 0 

in D, 

(3.4a) 

curl $ + ih: 2 ^ + X7(p = Uh 

in hi, 

(3.4b) 

div \I> = 0 

in D, 

(3.4c) 

$ x n - k\ I/ 4 = 0 

on dfl, 

(3-4d) 

ip = 0 

on dCl. 

(3.4e) 

Due to (p G Hq(Q), we easily obtain 

1M \i,n < C\\u h \\ 7h . 


(3.5) 

When the estimate (2.7) holds, taking J = Uh in (2.7) we have 


K ll^ll§+ Q ,n + ll curl ^lli+ Q ,n + «(1 + 

o,n < C(1 + k)\\ u h \ |j h , 

(3.6) 

which implies 

ll $ l|i+a,n < C(l + K)\\u h 

Ik- 

(3.7) 


By the equation (|2.6a[) with J = u hl we directly have 


(curl curl th, ^f) n - k 2 (^, i &)n + (iV<£>, ^)n = (i u h , ^)n, 
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which together with the fact (V</?, \1 /)q = 0 and the boundary condition (2.6c) yields 
(curl Tq curl ty) n - K 2 (Tq ^) n - (i/c 1 ^ 4 , = (i u h , ^)n- 


Thus, taking the imaginary part of the left-hand side of the above equation and using (3.6) 
we have 


«ll* llo,an s IKIkll'I'llo.n < Cn ||u k || Tk . 


(3.8) 


Next we present a key equality. 


Lemma 3.2. Let (<&, th, </?) be the solution of the dual problem We have 

6 


iKik=E r *’ 


k =1 


where 


T\ = (u{ x n - u h x n, $ - n v f)3T k , 

T 2 = « • n - u n h ■ n,ip — YlQ(p)o<j h , 

^3 = —{n(uh - u h)i ^ - nu*)d7 h , 

T a = -( Ku h + w h x n, T f< )aa, 

T 5 = -(f,Uu^) 7h , 

T 6 = (a h , div - n [ /T'))o' h - (a h , (T' - II v ^) ■ n) d7h . 


Proof. Using the dual first-order system (3.4), we obtain 


u h\\l h = (u h , curl $ + i/s; 2 ^ + Wif) 7h + (w h , i$ - curl ^) Th 

= (u h , curl n v $ + i/cn u^f + VIlQ^)j h + (w h , in v $ - curl 

+ (u h , curl ($ - n v $)K + (u h , i k 2 (^ - ny*))^ + (u h ,V((p - Ilg^))^ 

+ (ti7*,,i($ - n v f)) Tfc - (wh, curl (T' - rit/T'))^. (3.9) 


By the definitions of Ylu and Ily, we have (u/ l ,iK 2 (T r — Ilf/fl/))^ = 0 and (to/i, i($ — 
n^f))^ = 0. Integrating by parts and applying the property of the L 2 -projections yields 

(u h , curl ($ - n v $)) Th = (curl u h , $ - n v $) Th + (u h x n, $ - 

= (%xn,$- H v <&) d7h 

= (u{ x n, $ - n v $)ar h , (3.10) 


(u h ,v{(p - n Q <p)) 7h = —(div u h , cp — n Q <p) 7h + (v%-n,<p- n Q ip)m h 

= « -n,<p- n Q v)d7 h , (3.11) 


— ('Wh, curl (SH - IIc/^))^ = -(curl w h ,^F - Ilf/T')^ - (w h x n, ^ - 11^)^ 

= ~(Wh x W + (w* h x n, nt/^ajfc- (3.12) 


and 
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Taking r/, = Ily# in the equation (2.3a), noting that u h x n is continuous across each 
interior face and using the boundary condition (3.4d), we obtain 

Oft,, curl = (iw h ,n v $)j h + (u h x n,n v $)ar h 


= (iw h , riy$) Th + (u h X n. Ily# - $} 
= (iw h ,n v $) 7h + (u h x n, Ily# - $} 


97 h 
97h 


(u h , $ x n) 

(*&«**} 


an 


on- 


(3.13) 


Taking qh = I Iqip in (2.3c) and noting that u h ■ n is continuous across each interior face and 
p G Hq(Q) we have 


(“i., vn Q ip) 7t = (ui ■ n, n e v>)wi = (“ft - n, n Q ip - 


(3.14) 


We further take v h = nt/th in ( |2.3b[ ) to get 
- (wh, curl II[/T , )j /i 

= ~{w h x n Uu^g^ + (i/t 2 u h ,U u ^)j h - (cr h , div n t /T , ) T?i 
+ (d h . Ilt/T' • n) 97h - (/, Uu^) 7h 

= ~(wl x n - Tt(ul - u h ), n U^)gj h + (m 2 u h , Ilt/T')^ - (a h , div Ilt/T' - div v h)j, l 
+ {<7hi nt/th ■ n — 4/ • n)ffj h — (/, n t /T , ) a - h , (3.15) 

where the above second equality holds due to the fact that div V L = 0, dh is contin uous 
across each interior face and dh = 0 on £®. Inserting the above equalities (3.10)-(3.15) into 
the right-hand side of (3.9), we obtain the result. This completes the proof. □ 


We can now give the proof of Theorem 2.1 


Proof. (Proof of Theorem |2.1[ ) We derive the upper bounds for Tf, ■ • • , T 6 in Lemma 3.2 under 
the assumptions in Theorem |2.1[ By the Cauchy-Schwarz inequality, the approximation 
properties of standard L 2 -projections, the inequalities (3.5) and (3.6), we have 


-i,h. 


* fn 


Ti < C\\rf(u h - u h )\\grj h T t 5 (-)“||$|| i +Qi n < C\\rf(u h - u h )\\&j h T t 5 (-) Q (l + re)||«fc|hr fc , 

p 2 p 

I i h i 

T 2 = ~(r n (a h - d h ),p- II Q (p)&j h < C\\rf(a h - a h )\\d7 h Tn (-) 2 |M|i,n 

i i h i 

< C\\rf(a h - d h )\\ d7h Tf (-Y\\u h \\<j h , 

T 3 < C\\Tf(u{ - ^)|| aT?i r/(^) a ||^||i +Qin < C\\t?(u\ - u h )\\ a <j h Tf{^) a (l + l/K)\\u h \\ 7h . 


Taking r] h = Pm * hi ( T3e ) and using the property of the L-projection operator Pm on 
£f and the inequality (3.8) yields 


T,i = — (nu h +w h xn, P m^ )an = ~(g,PM)an < ||flf||o, aall^ ||o, an < Ck 
For the estimate of T 5 , we easily deduce 

T5 < H/llo.ull^llo.fi < C\\f\\ 0} nK 1 ||u/ l || a ' h . 


-1 


g\\o,dn\\Uh\\ 7h . 
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Applying integration by parts on T 6 , we have 

T& = — (Vcr/j, Sh — + (ah — <Jh, — n[/\E r ) ■ n)&x h = (ah — cr h , (fh 


• n) d7h 


< C\\tJ (a h - a h )\\d7 h Tn* (jell's?\\i +oli n < C\\T^(a h - 5\)IIot/T„ 5 (^)“(1 + 1/K)\\u h h h - 
Combining the above estimates for Ti, • ■ • , T 6 , we obtain 
llWill^y < Ck _1 (||/||o,n + llfl'llo.ao) llw^H^ 

+ c '( T t _5 (^)“( 1 + «) + r t (^) a (! + !/«)) lkt 5 («A - sDlKKIk 


+ C(tv ?(—)2 + r n 2 (-)"(l + l/n))\\r^(a h 

v p p ' 


i Y 


cr h)ll9T h ||^/i||^- 

Here we choose r t — j- and r n = (1+ ( ); ' . By the Young’s inequality, we have 


(3.16) 


\u h \\l h < c(k 2 


2 i K - 2 ||nll 2 

+ K l|y|lo,an 


0,0 


((! + 


2 A2a +1 + (l + l//c) 2 ( ^2 a -n 


K - 

P 


P 






+ ((i + K )(b 2 + h±d(-) 2 “- 1 )lldK 


^)ir 




Combining the above estimate and (3.1), the absolutely stable property of IIw^IIt^ can be 
easily observed, and we can further obtain (2.11) by the Young’s inequality in the regime 
nh/p > Cq. Then, by the estimate (3.2), we can also see the absolutely stable property of 
and have 

1 


\™h\\l h < 2ac 2 |K||^ 


+ Iloilo,n + “Iloilo,an- 


Then (2.12) is derived by (2.11). Furthermore, combining the fact that (cf. [25]) 

lKller h < Cph~^\\u h y h , 


(3.1), (2.11) and the triangular inequality yields the absolutely stable property of Ilff/JIdoy 
and the estimate (2.13). 

Whe n f — 0 and g — 0 in the first-order system (2.1), the estimates (2.11 )-(2.13) and 
Lemma 3.1 imply Wh — 0, Uh — 0 on 7h and u h = 0, oy = ah on d7h■ It then follows from 
(2.3b) that for any v h £ Uh, 

~(a h , div v h )j h + (a h , v h ■ n) d7h = (Va h , v h )j h - (a h , v h ■ n) d7h + (a h , v h ■ n) d7h = 0, 

which implies ah is piecewise constant on 7h . Due to the fact that ah — ah — 0 on dfl, we 
have ah = 0 on 7 h and cr^ = 0 on d7 h . Hence, the well-posedness of the HDG method (2.3) 
always holds without imposing any mesh constraint, i.e., the HDG method (2.3) is absolutely 
stable. □ 


Moreover, under the assumptions made in Theorem 2T, we can further get the upper 
bounds for ||curl Uh\\ j h and ||div Uh\\7 h - We take Vh = curl-tp, in (2.3a) to get 

(ittf h ,curlu h ) Tfc - (u h , curl curl u h )j h + (u h x n, curl u h ) a7h = 0. 


Using integration by parts on the above equation, we have 

(iw h , curl u h ) 7fi - (curl u h , curl u h )j h - (u[ x n, curl u h ) d7h + (u h x n, curlu^)^ = 0, 
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which directly yields 

1 _1 ± 

||curl < \\w h \\ 7 h \\cur\u h \\ 7h + C\\t? - u h )\\ d 7 h T t 2 ph-* ||curl u h \\ 7fi . 

Combining the above inequality, (3.1), (2.11) and ( |2.12[ ), then we get 

||curlMft||j h < c({- +p 2 + (1 + — )C s 2 tab ) ||/|| 0 ,n + (— + + (1 + — )C'stab)||fl r ||o,an)• 

\ K K K, 2 K, 2 h / 


Taking q h = divw ?t in (2.3c) and using integration by parts, we have 

1 i _i 

ul ■ n, div u h ) d7fi < C\\t£ (u h - a h )\\ 87 h rfph 2 ||div u k \\ 7h . 


|div u h \\ 7h = « 


■ n 


Then we obtain the upper bound for ||divit^||j h by the above estimate, (3.1) and (2.11) as 
follows: 

l|divui||j„ < C(1 + K)ipi ((1 + ^JH/Ilwi + (T + — )llsllo.ao' 

' K K. 2 


K 


Remark 3.1. By the estimates (3.1) and (2.11), we can get the upper bound for \\rf(u h — 
u h ) llaov Moreover, taking v h = Voy in (2.3b) and applying integration by parts, the Cauchy- 
Schwarz inequality, trace inequality and the estimates in Lemma 3.1 and Theorem 2.1, we can 


also get the stability estimate for 11 Voy 11 7h . When nh/p < Co, one may tune the parameters 
T t and r n according to the derivation of upper bound for the right-hand side of (3.16) and 
get the stability estimates. 


4. Error analysis 


In this section we provide detailed proofs of the a priori error estimates in Theorem 222 
We denote 


e w = Ylyw - 
e-u n = Pmu' 


w h , e~t = P M W* 


w 


h ? 


e u = Hull - u h , e~t = P M u l 


u 


h > 


U 


h > 


n Q <j oy, ej — l ) \/u Uh. 


In the following we first present the error equation for the analysis. 


Lemma 4.1. Let ( w,u,u ) and (wh,Uh,u h , cr h ,Uh) solve the equations (2.1) and (2.3). We 
have 


(i e w ,r h ) 7h - (e u , curl r h ) 7h + (e 2 t x n,r h ) d7h = 0, (4.1a) 

(e wj cur\v h ) 7h - (e-t x n,v h ) d7h + (m 2 e u , v h ) 7h 

- (e a ,d\vv h ) 7h + (e$,v h ■ n) d7h = 0, (4.1b) 

- (e u , Vq h ) 7h + (e G « • n, q h )&j h = 0, (4.1c) 

(e®*xn,ij h ) Wk \ 8 n = 0, (4. Id) 

(e®* x n, r) h )an + (/ce a «, 17 h ) d n = 0, (4.1e) 

( e u n ■ n, ^h)d 7 h = 0, (4. If) 


for all (■ r h ,v h ,r) h ,q h ,£ h ) eV h xU h x xQ h x M h ( 0). 


Proof. We notice that the exact solution (w , u , vf\ £ h , u, cr| £h ) also satisfies the equation (2.3). 
Hence, due to the property of standard L 2 -projection, the solutions Wh, Uh, u h \^ h , 

u'h\e. h , oy, uy|e,, in the equation (2.3) can be replaced by Ylyw, Pmw t \e h , LI U u, Pm^ 1 |g h , 
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PmU u |e h , U Q a,P M a\ eh respectively to derive a new equation, which subtracts the equation 
(2.3) to yield the result. □ 


Next we are going to present our first error estimate. 


Lemma 4.2. If we choose r t — f and r n = 


(1 -\-k)h 


we have 


1 

|«5e a t|| 0i en + ||r t 2 (e^ - e a t)|| aTh + ||r|(e CT - e 5 )|| aT ^ < Cr}(w,u), (4.2) 

lewllo-h < K||e u ||T h + C , 77('iu,u), (4.3) 


w/iere r)(w,u) = K\\w\\ t ,n + (l + (1 + n) 2)^||w|| s ,n, s > \,t > 


Proof. Let Vh = e w ,Vh = e u ,rj h = e^t,q h = e a ,fh = ej in the error equation (4.1). Then 
we get the following equalities after some simple manipulations which includes applying 
integration by parts: 


- (ie™, e w ) 7fi - (e w , curle u ) T?i + (e^ x to, <) OT(i - (e^ x to, e~t) d7h = 0, (4.4a) 

(e w , curle u )cr h - (e~t x TO,e^) OTh + (m 2 e u ,e u ) 7h - (e a ,div e u ) 7h + (e 5 ,e” ■ n)ffj h = 0, 

(4.4b) 

(e a , div e u ) 7h - (e CT , < • n) aTh + (e a , e^ ■ n) d7h = 0, (4.4c) 

(e& x n, e~t) d7hXm = 0, (4.4d) 

x n, e G t) an + («e a *, e s t) aQ = 0, (4.4e) 

- (e$, e G « ■ to) or* = 0. (4.4f) 


Adding the above equalities (4.4a)-(4.4f) together yields 


~(ie w ,e w ) 7h + ( in 2 e ui e u ) 7h + ((e^ - e~ t ) x n, < - e~t) d7h 

+ (Ke n t,e~t) dn - (e a - ea, « - ea") • n) 9Th = 0. (4.5) 


By the definition of Wh hr (2.4) we have 


( e w ~ e w t ) x n = (n v m - w h y x n - {P M w t - w h ) x n 

= (nvtu) 4 x n — x n — T t (u^ — u^) (4.6) 

= (n v it>) f x n - P M ™‘ x n - T t ((nc/W - e u ) 1 - ( P M u 1 - e a t)) 

= (n v wf x n - P M w t x n + r t (e^ - e G t) - r 4 ((n t/ w) f - P m m‘). 


Moreover, by the definition of uj) in (2.4), we have 


( e 


n 

it 


ea «)■n 


(LI uu — Uh) ■ n — ( P M u n - uf) ■ n 
(II uu - P M u n ) ■ n + r n (cr h - a h ) 

(II uu - P M u n ) ■ n + r n ((riQcr - e a ) - (P M cr - ea)) 
(II uu - P M u n ) ■ n ** T n (e a — e-s) + T n (n Q cr - P M cr). 


(4.7) 
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Inserting (4.6) and (4.7) into (4.5), we obtain 


+ m 2 \\e u f Jh + /c||e a t||g ian + || - e G t) f d7h + ||r ? f {e a - e 3 


2 + IIn an + II 


\&J h 


= -((n v iy ) 4 X n - P M w t x n, e* u e~t) d7h + {rffU-uu ) 1 - _ e u*)an h 

+ (e ff - e 5 , (II v u - P M u n ) ■ n) d7h + (e a - e d , T n (Ii Q a - P M o))m h 
= -((n v w -w)xn,e t u - e~t) d7h + (r t ( U v u - u), e* u - e~t) dTh 
+ (e a - e 5 , {Iluu - u ) ■ n)m h + (e CT - e?, r n (n Q o- - cr))^, (4.8) 

where the second equality is derived by the properties of L 2 -projections Pm and Pm- Based 


on (4.8), taking the real part and imaginary part of the left-hand side of (4.8) respectively, 


the estimates (4.2) and (4.3) can be obtained by the approximation properties of standard 
L 2 -projections, the Young’s inequality and the fact that cr = 0. This completes the proof. □ 

Now we start to use the duality argument to get an estimate for e u . Given e u e L 2 ( fl), 


we introduce the first-order system of the dual problem (2.6) with J = e 



i4> — curl 4/ = 0 

in fi, 

(4.9a) 


curl $ + ift 2 4/ + \7p> = e u 

in 12, 

(4.9b) 


div 4> = 0 

in 12, 

(4.9c) 


4> x n - = 0 

on dfl, 

(4.9d) 

Similar to the estimates in 

IMIi,o < C\ |eJ| Th , 

tp = 0 

(3.5)-(3.T) , we have 

on on. 

(4.9e) 

(4.10) 

ll^ll|+a,n + K \ ^ll§+a,0 

+ | curl th i +Q , i0 + (1 + k) 

k 4/| o,o < C( 1 + k)\ |e u | |j h . 

(4,11) 


Next we first present an important equality. 


Lemma 4.3. Let (4?,\E f ,(^) be the solution of the dual problem (f.9). It holds that 

5 


-u\\ 7h 


T, Ek-, 


k =1 


where 


Ei = (« - e G t) x n, $ - n v $>ar h , 

E 2 = (« - e G n) -n,ip- n Q<p)gn h , 

E-.i = ((e~t - elf) x n, - 11^$)^, 

Ea = -{e^t x n + /ce G t, 

P 5 = (e CT , div (\t r - n t7 T')) Th - (e 5 , (4/ - 11^) 




Proof. By the dual problem (4.9), we have 

IU II2 


(4.12) 


uiiT h = ( e «> curl $ + i/=c ^ + V<p) Th + (e w , i$ - curl ^) Tfi 

= (e u ,curin v $ + + VIIq^)^ + (e^illv^ - curin [7 T') Tft 

+ (e u , curl (4? - n v ^))T h + (e„,i/c 2 (^ - n [/ T')) Th + (e u , V(<p - IIq^))^ 

+ (e w ,i($ - n v $)) Th - (e w , curl (4/ - IV®')) Th . (4.13) 
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By the definitions of Yly and Ily, we have (e u ,iK 2 (\E r — Lljy'I'))^ = 0 and (e w ,i($ — 
n v $)) T;> = 0. Similar to the derivations of (3.10)-(3.12), we have 

(e u , curl ($ - n v $)) T/i = (el x n, $ - IlyS)^, (4.14) 

(e«, V(<P ~ n q<p))t h = (e” • n,ip - H Q p) d7h , (4.15) 

~(e w , curl (^ - Uu^)) 7h = -(el x n, ^ - n^)^. (4.16) 


Taking rg = Ily# in the equation (4.1a), noting that e a t x n is continuous across each 
interior face and using the boundary condition (4.9d), we obtain 

(e u ,curin v $K = (ie^.IIv^K + (e fi t x n.Il v $) d7h 

= (ie„„ n v $) Th + (e G t x n, 11^4? — *&)dr h — (e^t, $ x n) 90 


— (ie«,, Ily#)^ + (e G t x n, IIy4? — &)&i h ~ (^e*, 


an- 


(4.17) 


Note that • n is continuous across each interior face and p G Hq(D,). We let qh = f IqP 
in (4.1c) to obtain 

(e u , Vn Q p) 7h = (e G " . n , H Q p) d7h = ( e ^ ■ n, Ii Q p - p) d7h . (4.18) 


We further take v h = rif/th in (4.1b) to get 

- (e w , curin t /T f )j, i 

= —( e w { x n, Hu^)g 7h + ( in 2 e u , II[/^ , )‘j h — (e<r, div Ilf/4/)^ + (eg, rit/T' • n)g 7h , 

= —(e-w* x n ? nr/T' — ^)d7 h ~ (e^t x n, 4 > 4 )qq + (ift 2 e u , Hu^) 7h 
- (e a , div II [/ T' - div 4>) Th + (eg, hlf/T' ■ n — 4> ■ n) d7h , (4.19) 

where the above second equality holds due to the fact that div 4/ = 0, e~ t x n and e g are 
continuous across each interior face, and eg = 0 on Ef. Then, inserting (4.14)-(4.19) into 
(4.13) yields the result. □ 

Based on the above lemma, we can obtain the estimate for ||e u ||T h . 


Lemma 4.4. If the regularity property ( f.10) holds, and r t ,T n are chosen as in Lemma 4-2, 
we have 


-uWVh 


< C(R 


wUI^II l+a,n 


Ru\\u\\l + a } n 


). 


(4.20) 


where R w and R u are defined as in Theorem 2.2 


Proof. We need to derive the upper bounds for E\, ■ ■ ■ , E§ in Lemma T3 By the Cauchy- 
Schwarz inequality and the approximation property of Ily, we obtain 

Ei < C||d(el - e a 0|| OTl rry-)l#||.«, n < C\\r}(e‘„ - e a .)|| OT ,(l + K )(-) l+ “IMk- 

P P 


By the identity (4.7) for (e™ — eg") • n and the fact that a = 0, we can derive that 
E 2 = ((n v u - u + u n - P M u n ) ■ n - T n (e a - eg) + r n (n Q a - P M cr), p - U Q p) d7h , 


<c(f) 

P 
Ji 


,h. 


i,n 


«l|i+a,n + Tj\\Tf(e a - eg)\\ d 7 h )(-fi\ 

p . :: 3 + (1 + p k); Y Ik - 


<G((-)i + “|| U ||. 
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Moreover, by the identity (4.6) for (e^ — e-t) x n and the triangular inequality, we get 
E 3 = —{i't{e t u - e a «), ^ - n r/ ^)aT, + {Pmw 1 x ra - (n^w)* x n, ^ 

+ ^.{(JluuY — Pm u f , v h — 


< C ( { - l±dA«-l||d(el-ei 

rC J) 


. (1 + k) . h^2 a \\ II 

U + IMIi +a ,n 

p 2 


K 


(1 + k) h s 


2a- 1 1 


U i 


K 


P 


2 




By the boundary condition (4.1e), we have = — (e~t x n + Ke a i, -Pm^ }ao = 0. Applying 


integration by parts, we obtain the estimate for E 5 as follows: 

E 5 = (Ve ff , ^ - n^) + (e a - e 5 , • n) d7h 

= {e a - e ? , (* - 11^) • n)ffj h < C { - 

hi p 

Finally, combining the above estimates for A], • • • , E 3 and the estimate (4.2), we can conclude 


)“" 5 || rl{e„-e- n 


\d7 h \\e u \\7 h - 


the result. 


□ 


We can now give the proof of Theorem 2.2 


Proof. (Proof of Theorem 2.2) By the triangular inequality, we have 

11“ - “fclk < 11“ - U uu Ik + ||e„|k- 


The error estimate (2.14) can be obtained by the approximation property of Tljj and the 


estimate (4.20) for ||e u ||o- h . Similarly, (2.15) can be obtained by the triangular inequality, 
the approximation property of IIv, (4.3) and (4.20). □ 


Remark 4.1. Besides we get the error estimate for ||(e CT — e$)\\gy h in (4.2), we can also 


obtain the error estimate for 11Ve o-Jlyry - Actually, this can be similarly derived as the stability 
estimate for || (cf. Remark [371 ) by taking Vy = Ve a in the error equation (4.1b). Then 

the error estimate for ||V(er — cr^jfk can be further deduced by the triangular inequality. 
When nh/p < Co, one may tune the parameters r* and r n (cf. Remark 3.1) and get the error 
estimates. 


5. Stability and error estimates for ideal case 


In this section, we consider the stability estimates and error estimates of the HDG method 
(2.3) under some ideal assumptions of the problem (1.1) and the dual problem (2.6). We 


assu me th at when D is a smooth star-shaped domain, the solutions of the first-order sys¬ 
tem (2.1) satisfy that u G H 2 (Q) and w G iT 1 (D). When / is divergence-free and 

— - " i 

g G Hf(dfl) := {g G [FC (<9D)] 3 , g ■ n — 0 on <9D}, we assume the following estimate 
holds, which has been mentioned in nm that 


IMkn + IH|i,n < C(1 + K)M{f,g) + C\\g\\i dQ , 


where M(f,g) = \\f\ 


the dual problem (2.6 
there also holds 


o,n + H^llo,an- In this ideal case, we can also assume the solution of 
satisfies that 4/ G iT 2 (D), the estimate (2.7) holds with a = | and 


ll^lkn < C{ 1 + «)|| J - VpHo.n < C(1 + k)|| J|| 0 ,n. 


(5.1) 
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We assume the approximation results of L-projections in (2.2a)-(2.2a) still hold, then we 
have the following stability estimates and error estimates for the HDG method (2.3). 


Lemma 5.1. We assume that (2.1) holds with a — \ and (5.1) also holds true. Let 
(wht Ufi, tt^, Chi &h) be the solution of the problem (2.3). We have 


k^lk < II — ||on + C s tab|| — ||o,an 

V K K 

\ W hh H < C((~ + G s 2 tab )||/|| 0 ,O + (— + C'stab) ||fl r ||o,Sf 2 ) , 

\ 1 % R2 

^s 2 tabll“llo,0 + C'stab||-||0,90 I, 


, „ /, nh ,i , i 

\u h \\dT h < c ((y ) 5 +P h 2 


K 


K 


(5.2) 

(5.3) 

(5.4) 


where := 1 + & 

Proof. In order to get the upper bou nd fo r Ikfel k, in deed, it only needs to b oun d the terms 

When (2.7) holds with a = \ and (5.1) also holds, 


2.1 


Tj, • • ■ , Tq as in the proof of Theorem 

we have the following regularity estimate for the dual problem (3.4), 

||3||i,n + ||^II 2,0 + kl^lkn + Ilcurl^H^n + k( 1 + «)||\I'||o,n < C( 1 + «)|kfclk- 
By the above regularity estimate, we have 

Ti < C\\rf (u l h - Wh)||ar h r t _5 (^)3(l + kIKIk; 

T 3 < C'||r/« - «fc)||ar h T t 5 (^)i(l + kIMk, 

T 6 < C\\Tn (<?h - O h )\\d7 h Tn 2 (“)^ (1 + K )\\ U hh h 

and the estimates for T^T^Ts are the same as the estimates in the proof of Theorem 2T 
Combining the estimates for Tf, • ■ ■ , T 6 again, we obtain 

IKIk < C'k ' 1 (H/||o,n + llflfllo,an) Ikfclk 

+ C'(r t _5 (^ +r/(^)i)( 1 + k)||t/« - u h )\\d7 h \\u h \\ 7h 

+ C(tu + r n 2 (^)^(1 + «))|| Tn{a h - k) ||ar fe || u h \\ 7h . 


(5.5) 


Here we choose r t — f and r n = 


_ (1 -\-K,)h 


Then, the stability estimate (5.2) can be obtained 


h ' p 

by (5.5), (3.1) and the Young’s inequality, and the stability estimates (5.3) and (5.4) can be 
further derived as the analysis in Theorem 2.1 □ 


Lemma 5.2. We assume that (2.1) holds with a — \ and (5.1) also holds true. Let 
(■ w h ,u h ,u h ,(Jh,(?h ) be the solution of the problem (2.3). We have 

Ik Y C* ((Riv ||||T Ru |k||s,n) j 


\w - w h Ik < C[ (— + kR 


s-1 


\W 


i h 

|t,a + ((1 + (1 + k 5 ) i + kR-v 

p s ~ 1 


where s > 1, t > 1, R w := (1+ ^+i fcf+1 H- ~t + T 


(l+Ac)/l t + 1 

pt 


and R u := ^ + 


( 1 +K)?h s 


+ 


klU.n 

(I+k)/*' 


(5.6) 

(5.7) 
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Proof. When (2.7) holds with a = \ and (5.1) also holds, we have the following regularity 
estimate for the dual problem (4.9), 


H^llpn + ||^II 2 ,o + + ||curlfE'IIgn + (1 + ft)ft||tE , ||o,f 2 < 0(1 + ft)||e u || 7 h . 

Repeating the similar estimates in Lemma |4.4 for Ei, ■ ■ ■ , E 5 and using (4.2), we obtain 

||e u ||T fc — O (f^w 11^l|t,n || u ||s,o)• (5.8) 


Then the error estimate (5.6) is obtained directly by the triangular inequality, the approxi¬ 


mation property of Ylu and the above estimate. The error estimate (5.7) can also be obtained 
by the triangular inequality, the approximation property of 11^, (4.3) and (5.8). □ 


Remark 5.1. By Lemma 5.2, under the assumptions made in the section, we have 

,^2^2 _ 


\u 


II <" r ( Kh2 
u hh h < C{ — 

\w-w h \\ Th < C(~~ ~ 


hi 


pz 

2 h 2 




2 h 2 


pz 


P‘ 

n^h 

p2 


jr 


P 


), 


+ —)M(f,g < C(— + — . 

pi / p pi 


Here n > 1 and M (/, g ) = M ( f.g ) + 1||gr ||i dCl . The above estimates indicate that the error 
11 w — Wh\\j h can not be controlled by —, and the pollution term is of order O(*-?-). This 
provides evidence of the existence of the so-called “pollution effect”. When WL < C, the 
discrete stability estimates for ||M h ||j h and || , iu/ 1 ||t, j can be improved as Hu/JI^ < — M(f,g) 
and ||io fc || Th < CM(f,g). 


6. Numerical results 

In this section, we present numerical results of the HDG method for the following time- 
harmonic Maxwell problem (cf. [TO] ) in a unit cube = [0,1] x [0,1] x [0,1]: 

curl curl u — n 2 u = 0 


curl u x n — inu* = g 


in D, 
on dfl. 


Here g is chosen such that the exact solution is given by 


— (J KZ J KX JKy\T 


u — e 


r- 


The time-harmonic Maxwell problem (1.1) is an approximation of electromagnetic scatter¬ 


ing problem with time dependence e lujt , where cu i s freq uency. If the problem is proposed with 
time dependence e~ IUJt , then the sign before i in (1.1b) is negative. The analysis of the HDG 


method in this paper fits well for both of cases. In the following experiment, we apply the 
HDG method with piecewise linear (HDG-P1), piecewise quadratic (HDG-P2) and piecewise 
cubic (HDG-P3) finite element spaces respectively to the second case. For the fixed wave 
number hi, we first show the dependence of the convergence of ||m — tt/i||o,n an d ||m — rCfe||o,o 
on polynomial order p and mesh size h. Figure [T| displays the above errors for n = 20 by 
the HDG-P1, HDG-P2, and HDG-P3 approximations. The pollution errors always appear 
on the coarse meshes. However, we find that the errors converge almost in 0(nh 2 /p 2 ) on the 
fine meshes, which is a little better than the theoretical prediction for ||iu — in/ t ||o,n- On the 
other hand, for the cases of n = 30 and n = 50, Figure [2] shows that the errors of ||w — W/i||o,n 
and ||in — iu/i||o,n always decrease for high order polynomial approximations. 
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HDG—PI (k=20) 



HDG—P2(k=20) 



HDG—P3(k=20) 



FIGURE 1. Errors of ||k, — Mft,||o,n and ||u> — xt?/,,11o.fi for k = 20 by the HDG-P1, HDG-P2 
and HDG-P3 approximations. 






FIGURE 2. Errors of ||m — iifc||o,fi and ||w; — tw/dlo.n for k = 30 and k = 50 by the 
HDG-P1,HDG-P2 and HDG-P3 approximations. 


Figure [3] displays the relative errors ||w — U/i||o,n/IM|o,n and ||iw — ||o,n/Iloilo,n for the 
HDG-P1 approximation according to different mesh size conditions. The left graph of Figure 
[3] shows the relationship between the relative errors and the wave number k under the mesh 
condition kIi = 2 for the HDG-P1 approximation. We observe that the relative errors cannot 
be controlled by k h and increase with k, which indicates the existence of the pollution error. 
The right graph of Figure [3] shows the relative errors of the fiDG-Pl approximation under 























































HDG method for Maxwell equations 


19 


the mesh condition tc’h 2 = 2. It shows that under this mesh condition, the relative errors 
do not increase with k. 




Figure 3. Left: The relative errors ||tt - tt/i||o,n/IM|o,n and \\w - Wh||o,n/|M|o,n for 
the HDG-P1 approximation under the mesh condition nh = 2. Right: The relative errors 
11“ - 'U/t||o,n/|Mlo,n and ||u; - tt7/j||o,n/|Mlo,n for the HDG-P1 approximation under the 
mesh condition n 3 h 2 = 2. 


For fixed wave number k, we show the relative error ||to — iUh||o,n/||'R ? ||o,n for the HDG- 
P1 approximation with respect to the relative error ||curl(it — ii/j)||o,n/||curlii||o,n for the 
standard lowest-order edge element approximation of the second type. The left graph of 
Figure [4] displays the relative error of the HDG-P1 solution for k = 10,20,30, while the 
right one shows the relative error for the same cases based on the standard lowest-order 
edge element method. We find that the relative error for the HDG-P1 approximation stays 
around 100% while the relative error for the standard edge element approximation oscillates 
around 100% before they are less than 100%, which confirms the stability property of our 
theoretical analysis for the HDG method and indicates that the HDG method is more stable 
than the standard edge element method for the time-harmonic Maxwell problem with high 
wave number. 



FIGURE 4. The relative error ||u; — iu/ t ||o,n/||'u>||o,n (left) for the HDG-P1 approxima¬ 
tion and the relative error ||curl(n — wDI|o,n/||curlu||o,n (right) for the lowest-order edge 
element (the second type) approximation for k = 10, 20,30 respectively. 

Table[l]shows the numbers of degrees of freedom (DOFs) and the relative error || curl (it — 
w/ l )||o,o/||curlu||o,n for the standard edge element approximation with respect to the relative 
error ||iy—iUh||o,n/||' iy ||o,o f° r th e HDG-P1, HDG-P2 and HDG-P3 approximations. It can be 
observed that the HDG-P1 approximation performs better than the standard edge element 
method when the numbers of DOFs are close. We can also find that the HDG method 
with higher order polynomial approximation may reach more accurate solutions with less 
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TABLE 1. The relative error ||curl(it — Wfc,)||o,fi/|| cur lti||o,n for the lowest-order edge 
element (the second type) approximation for the case k = 50 and the relative error ||io — 
^/i||o,n/||i«||o,n for the HDG-P1, HDG-P2 and HDG-P3 approximations with respect to 
different DOFs. 



DOFs 

8368 

62048 

477376 

3744128 

Edge element 

The relative error 

111.9% 

115.8% 

109.1% 

42.7% 


DOFs 

9792 

76032 

599040 

4755456 

HDG-P1 

The relative error 

96.8% 

96.7% 

82.3% 

30% 


DOFs 

— 

19584 

152064 

1198080 

HDG-P2 

The relative error 


100% 

89.4% 

21.6% 


DOFs 

— 

32640 

253440 

1996800 

HDG-P3 

The relative error 

— 

100% 

54.3% 

2% 


DOFs, which indicates the efficiency of the HDG method with high polynomial order for 
the time-harmonic Maxwell problem with high wave number. We should note that the 
numerical results in [10] show the stability of the IPDG method based on the piecewise 
linear polynomial approximation for the time-harmonic Maxwell problem with high wave 
number. Here, our HDG method preserves the advantages of the IPDG method in [10], and 
it results in a discrete system with significantly reduced DOFs when it is applied for the high 
order polynomial approximation. 





FIGURE 5. The traces of the real part of the first component of the HDG-P1 solutions 
for k = 30 and k = 50 (left and right) on the meshes with h = 1/16 and h = 1/32 (top and 
bottom). The traces of the real part of the first component of the exact solution are plotted 
in the blue lines. 

For more detailed comparison between the HDG methods with different polynomial order 
approximations, we consider the problems with wave number k = 30, 50. We restrict the 
solution plot in the line segment {(x,y,z) : x = 0.5, y = 0.5,0 < z < 1} and observe the 
traces of the real part of the first component of the HDG solutions. The traces of the real part 
of the first component of the exact solution are also plotted in the blue lines in Figure [5] and 
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FIGURE 6. The traces of the real part of the first component of the HDG-P2 and HDG- 
P3 solutions (top and bottom) for re = 30 and re = 50 (left and right) on the mesh with 
h = 1/16. The traces of the real part of the first component of the exact solution are plotted 
in the blue lines. 


Figure [6] The left graphs of Figure [5] display the traces of the real part of the first component 
of the HDG-P1 solution on the meshes with h = 1/16 and h = 1/32 for k = 30, while the 
right graphs of Figure [5] show the same traces for k = 50. Figure [6] displays the traces of the 
real part of the first component of the HDG-P2 and HDG-P3 solutions on the mesh with 
h = 1/16 for k = 30,50 (left, right). On the coarse mesh with h = 1/16, the shapes of the 
HDG-P2 and HDG-P3 solutions are roughly the same as the exact solution while the shape 
of the HDG-P1 solution does not match the exact solution well. We can also observe that the 
HDG solutions of high order polynomial approximations on the mesh with h = 1/16 perform 
even better than the HDG-P1 solution on the mesh with h — 1/32 especially for k = 50, 
which shows the advantage of the HDG method with high order polynomial approximation 
for the time-harmonic Maxwell problem with high wave number. Thus, although the phase 
error appears in the cases of coarse mesh and low order polynomial approximation, it can 
be reduced in the fine meshes or by high order polynomial approximations. 
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